library(DESeq2)
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(pheatmap)
library(outliers)
library(dplyr)
library(plotly)
library(grid)
library(futile.logger)
library(VennDiagram)
options(width=120)

В ходе данной домашней работы понадобятся те же файлы, что и в лекции: “GSE89225_illumina_counts.csv”, “conditions.csv”, “human_mart.txt”. Для начала убедимся в том, что мы можем эти файлы прочитать. И посмотрим, что в них находится.

Неуместного аутлайера с условным обозначением treg_NBP_patient3 я вычислила, но немного позднее по ходу этого файла. Но, хотя Читатель узнает его имя позднее, я, раз знаю немного больше, удалю его из данных прямо сейчас сразу после считывания файла.

counts <- read.csv("GSE89225_Illumina_counts.csv", row.names=1)
conditions <- read.csv("conditions.csv", row.names=1)
mart <- read.table("human_mart.txt", sep="\t", header=1, check.names = F)

# а вот и его ненастоящее имя treg_NBP_patient3

counts <- select(counts, -treg_NBP_patient3)

conditions <- conditions[-12,]

# а это несколько строчек, чтобы взглянуть на струкутру датафреймов 
print(counts[1:6, 1:2])
##                 tconv_tumor_patient10 cd177negtreg_tumor_patient10
## ENSG00000000003                   510                          935
## ENSG00000000005                    25                            1
## ENSG00000000419                   613                          818
## ENSG00000000457                   269                          587
## ENSG00000000460                    55                          142
## ENSG00000000938                    36                           69
dim(counts)
## [1] 63677    33
head(conditions)
##                                            tissue                               cells
## tconv_tumor_patient10        tissue: breast tumor cell type: Conventional CD4 T cells
## cd177negtreg_tumor_patient10 tissue: breast tumor       cell type: Regulatory T cells
## cd177postreg_tumor_patient10 tissue: breast tumor       cell type: Regulatory T cells
## tconv_tumor_patient1         tissue: breast tumor cell type: Conventional CD4 T cells
## treg_tumor_patient1          tissue: breast tumor       cell type: Regulatory T cells
## tconv_NBP_patient2                    tissue: NBP cell type: Conventional CD4 T cells
dim(conditions)
## [1] 33  2
head(mart)
##           Gene ID Associated Gene Name Gene type EntrezGene ID
## 1 ENSG00000252303            RNU6-280P     snRNA            NA
## 2 ENSG00000281771                Y_RNA  misc_RNA            NA
## 3 ENSG00000281256         RP11-222G7.2   lincRNA            NA
## 4 ENSG00000283272      Clostridiales-1      sRNA            NA
## 5 ENSG00000280864        RP11-654C22.2   lincRNA            NA
## 6 ENSG00000280792        RP11-315F22.1   lincRNA            NA
dim(mart)
## [1] 64101     4

RNA-seq

Немного слов учителя об РНК-секвенировании:

Rna-seq steps:

  • Изоляция РНК
  • Rna selection / depletion
  • Фрагментация
  • Синтез кДНК
  • Секвенирование

Rna selection / depletion:

  • вся РНК
  • тянем за поли(А)-хвосты (только мРНК)
  • удаляем рибосомальную РНК (ribo-zero kit)
  • таргетное секвенирование

Why Rna-seq?

  • Не ограничены существующей сборкой и дизайном микрочипа
  • Низкий фоновый сигнал
  • Точность позволяет смотреть экспрессию отдельных изоформ

Sanity checks

###…и необходимой чистоте данных

Нужно всегда проверять длины библиотек и количество rRNA reads, которые оказались в библиотеке. Количество ридов можно проверять после выравнивания или после квантификации.

# Повторим сказанное на языке машины

proteinCoding <- mart[mart[, 3] == "protein_coding", ]
rRNA <- mart[mart[, 3] == "rRNA", ]

pcCounts <- counts[rownames(counts) %in% as.character(proteinCoding[, 1]), ]
rrnaCounts <- counts[rownames(counts) %in% as.character(rRNA[, 1]), ]

sampleCount <- ncol(counts)
toPlot <- data.frame(
  sample=rep(colnames(counts), 3),
  value=c(colSums(counts) - colSums(pcCounts) - colSums(rrnaCounts), 
          colSums(pcCounts), 
          colSums(rrnaCounts)),
  type=c(rep("other", sampleCount), 
         rep("protein coding", sampleCount),
         rep("rrna", sampleCount))
)

plot <- ggplot(data=toPlot, aes(x=sample, y=value, fill=type)) +
  geom_bar(stat="identity") + theme_bw() + 
  theme(axis.text.x = element_text(angle=90, vjust=0.5))
plot

DESeq2

DESeq2 и правда хорош и, возможно, даже красив. Тут и дифференциальная экспрессия, и нормализации, и PCA-plots.

# далее я создаю целых три DESeq-объекта, мне показалось, что так удобнее (вообще-то сначала нет, я пыталась работать с аргументом contrasts). Права ли я?

dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = conditions,
                              design = ~ tissue + cells)

dds <- dds[rowSums(counts(dds)) > 20, ]
dds <- DESeq(dds)
vst_dds <- vst(dds)
counts.norm <- assay(vst_dds)

# а вот и тот самый вышеупомянутый outlier, встречаем: 

outlier(counts.norm)
##        tconv_tumor_patient10 cd177negtreg_tumor_patient10 cd177postreg_tumor_patient10         tconv_tumor_patient1 
##                     18.31254                     19.06682                     18.75782                     19.50822 
##          treg_tumor_patient1           tconv_NBP_patient2            treg_NBP_patient2         tconv_tumor_patient2 
##                     19.42026                     18.74975                     18.97194                     18.68064 
##  cd177negtreg_tumor_patient2  cd177postreg_tumor_patient2           tconv_NBP_patient3         tconv_tumor_patient3 
##                     18.86007                     19.00433                     18.76385                     19.09414 
##  cd177negtreg_tumor_patient3  cd177postreg_tumor_patient3           tconv_NBP_patient4            treg_NBP_patient4 
##                     18.73277                     19.78128                     18.74858                     19.24072 
##         tconv_tumor_patient4  cd177negtreg_tumor_patient4  cd177postreg_tumor_patient4          treg_tumor_patient5 
##                     18.51430                     19.73984                     19.74446                     18.87674 
##         tconv_tumor_patient5         tconv_tumor_patient6  cd177negtreg_tumor_patient6  cd177postreg_tumor_patient6 
##                     18.84609                     18.85645                     18.68582                     18.93730 
##           tconv_NBP_patient7            treg_NBP_patient7           tconv_NBP_patient8            treg_NBP_patient8 
##                     18.55074                     18.03944                     18.73227                     18.73946 
##          treg_tumor_patient8         tconv_tumor_patient8           tconv_NBP_patient9            treg_NBP_patient9 
##                     18.43248                     18.40272                     18.90248                     18.56628 
##         tconv_tumor_patient9 
##                     18.94650
max(outlier(counts.norm)) # а это treg_NBP_patient3
## [1] 19.78128
# Вот такой еще

dds_tissue <- DESeqDataSetFromMatrix(countData = counts,
                              colData = conditions,
                              design = ~ tissue)

dds_tissue <- dds_tissue[rowSums(counts(dds_tissue)) > 20, ]
dds_tissue <- DESeq(dds_tissue)
vst_dds_tissue <- vst(dds_tissue)
counts.norm_tissue <- assay(vst_dds_tissue)

# И вот такой еще

dds_cells <- DESeqDataSetFromMatrix(countData = counts,
                              colData = conditions,
                              design = ~ cells)

dds_cells <- dds_cells[rowSums(counts(dds_cells)) > 20, ]
dds_cells <- DESeq(dds_cells)
vst_dds_cells <- vst(dds_cells)
counts.norm_cells <- assay(vst_dds_cells)

Поглядим на PCA без выброса?

pca_data_tissue <- prcomp(t(counts.norm_tissue))
percents_tissue <- pca_data_tissue$sdev^2 / sum(pca_data_tissue$sdev^2)
to_plot_tissue <- t(counts.norm_tissue) %*% pca_data_tissue$rotation

pca_data_cells <- prcomp(t(counts.norm_cells))
percents_cells <- pca_data_cells$sdev^2 / sum(pca_data_cells$sdev^2)
to_plot_cells <- t(counts.norm_cells) %*% pca_data_cells$rotation

gdata_tissue <- data.frame(
  x=to_plot_tissue[, 1],
  y=to_plot_tissue[, 2],
  tissue=conditions[, 1],
  name=rownames(conditions)
)

gdata_cells <- data.frame(
  x=to_plot_cells[, 1],
  y=to_plot_cells[, 2],
  cells=conditions[, 2],
  name=rownames(conditions)
)

# немного ненавязчивого интерактива от plotly

ggplotly(ggplot(data=gdata_tissue, aes(x=x, y=y, color=tissue, shape=tissue, text=name)) +
  geom_point(size=3) + theme_bw()  +
  xlab(paste0("PC", 1, ": ", formatC(100 * percents_tissue[1], digits=4), "%")) +
  ylab(paste0("PC", 2, ": ", formatC(100 * percents_tissue[2], digits=4), "%")))
ggplotly(ggplot(data=gdata_cells, aes(x=x, y=y, color=cells, shape=cells, text=name)) +
  geom_point(size=3) + theme_bw()  +
  xlab(paste0("PC", 1, ": ", formatC(100 * percents_cells[1], digits=4), "%")) +
  ylab(paste0("PC", 2, ": ", formatC(100 * percents_cells[2], digits=4), "%")))
ggplotly(plotPCA(vst_dds, intgroup=c("tissue", "cells"), ntop = 500) + theme_bw())

Differential expression

res <- results(dds)
res <- res[order(res[,4]),]

res_tissue <- results(dds_tissue)
res_cells <- results(dds_cells)

res_tissue <- res_tissue[order(res_tissue[, 4]), ]
res_cells <- res_cells[order(res_cells[, 4]), ]

“Самый обычный способ визуализировать данные дифф.экспрессии – это volcano plot. По оси x мы будем откладывать log fold change, а по y - adjusted p value”.

Попробую обличить слова в действительность.

genes_tissue <- rownames(res_tissue)
genes_cells <- rownames(res_cells)

gdata_tissue <- data.frame(
  log_fold_change=res_tissue$log2FoldChange,
  p_adj=res_tissue$padj,
  n = "Breast tumor vs Normal breast tissue"
)

gdata_cells <- data.frame(
  log_fold_change=res_cells$log2FoldChange,
  p_adj=res_cells$padj,
  n = "Treg vs Tconv"
)

# я все это делаю еще и потому, что потом мне отдельно понадобятся датафреймы для клеток и тканей

# но для volcano, хотя я и планировала сначала использовать grid_arrange (почему и создавала два gdata для клеток и для тканей), в конечном итоге я использую слитый воедино датафрейм general_gdata. Потому как в последующем применю facet_grid.

general_gdata <- rbind(gdata_tissue, gdata_cells)
general_res <- rbind(res_tissue, res_cells)

# создадим порог значимости

general_gdata$threshold <- general_gdata$p_adj <= 0.01
general_gdata$sig <- sapply(general_gdata$threshold, function(x) ifelse(x == TRUE, "Significant", "Not Significant"))
general_gdata <- na.omit(general_gdata)

# Вы только посмотрите, какая красота!

ggplot(general_gdata,aes(x=log_fold_change, y=-log10(p_adj), col = sig)) +
       facet_grid(. ~ n) +
       geom_point(size=1) +
       theme_bw() +
       scale_colour_manual(values = c("black", "red")) + 
       labs(x = "Log fold change", y ="Adjusted p.value")+
       geom_hline(yintercept=2, linetype = 2, size = 1, colour = 'red')

Поскольку я задержалась с домашкой, я решила, а может, еще что-нибудь почитать и как-нибудь поиграть с построением volcano? Прием, которым я воспользовалась ниже, направлен на более тонкое выявление дифференциальной экспрессии, основанное и на значении log2 fold change между двумя уровнями факторной переменной (уровнем экспрессии генов в разных клетках и в разных тканях) и на значении p-adj (p-val после поправки BH (Benjamini & Hochberg (1995)). Из полученного графика мы видим, что хотя различия в экспресии по показателю log2 fold change между некоторыми генами и выявляются, они остаются тем не менее незначимыми (черным цветом).

p_adj <- c(general_res$padj)

log_fold_change <- c(general_res$log2FoldChange)

this_frame <- general_gdata[,-c(4,5)]

this_frame_green <- subset(this_frame, log_fold_change < -1 & p_adj < .01) # определим зеленую зону

this_frame_green <- cbind(this_frame_green, rep(1, nrow(this_frame_green)))

colnames(this_frame_green)[4] <- "Color"

this_frame_black <- subset(this_frame, (log_fold_change >= -1 & log_fold_change <= 1) | p_adj >= .01) # определим плохую черную зону

this_frame_black <- cbind(this_frame_black, rep(2, nrow(this_frame_black)))

colnames(this_frame_black)[4] <- "Color"

this_frame_red <- subset(this_frame, log_fold_change > 1 & p_adj < .01) # определим красную зону

this_frame_red <- cbind(this_frame_red, rep(3, nrow(this_frame_red)))

colnames(this_frame_red)[4] <- "Color"

frame_total <- rbind(this_frame_green, this_frame_black, this_frame_red)

frame_total$Color <- as.factor(frame_total$Color)

genes <- rownames(general_res)

##Сконструируем объект volcano

my_another_volcano_hello <- ggplot(data = frame_total, aes(x = log_fold_change, y = -log10(p_adj),col= Color)) +
  geom_point(alpha = 0.5, size = 1) + theme_bw() + theme(legend.position = "none") +
  geom_hline(yintercept = 2, colour = "black", linetype = 2, size = 1) + 
  scale_color_manual(values = c("green", "black", "red")) +
  facet_grid(. ~ n) +
  labs(y="Adjusted p.value", x="Log fold change")

ggplotly(my_another_volcano_hello)

А теперь взглянем на горячую карту. Горячая карта покажет горячие точки. Вдруг эти точки о чем-то говорят. (Шучу, - конечно же, они говорят)

# процедуры извлечения генов из пасвея

kkeys <- keys(org.Hs.eg.db, keytype="ENSEMBL")
goAnno <- AnnotationDbi::select(org.Hs.eg.db, keys=kkeys, 
                                keytype="ENSEMBL", columns=c("GOALL", "ONTOLOGYALL", "SYMBOL"))
goAnno <- tbl_df(goAnno)
goAnno <- filter(goAnno, GOALL=="GO:0007159")
# or you can pick ENTREZ, or SYMBOL, or whatever you want
genesToVisualise <- goAnno$ENSEMBL

# несколько шагов перед тем, как предложить что-то изящной функции

to_visualise_zzzz <- counts.norm[rownames(res), order(conditions[, 2])]

subsetting <- rownames(subset(res, rownames(res) %in% genesToVisualise))

nearly_want_to_see <- to_visualise_zzzz[subsetting,]

exactly_want_to_see <- t(apply(nearly_want_to_see, 1, function(r) {(r - min(r)) / (max(r) - min(r))}))

# О, как она изящна

pheatmap(exactly_want_to_see, 
         show_rownames = F, cluster_rows = F,
         cluster_cols=F,
         annotation_col = conditions,
         main = "GO:0007159: leukocyte cell-cell adhesion")

А вот здесь я получаю пересечение размером 1. Burn the witch?

give_me_tissue <- subset(general_gdata, n == "Breast tumor vs Normal breast tissue")
give_me_cells <- subset(general_gdata, n == "Treg vs Tconv")

what_want_from_tissue <- sapply(give_me_tissue$p_adj, function(x) ifelse(x < 0.01, give_me_tissue$sig, "No"))

what_want_from_cells <- sapply(give_me_cells$p_adj, function(x) ifelse(x < 0.01, give_me_cells$sig, "No"))

finally_desire_tis <- what_want_from_tissue[what_want_from_tissue == "Significant"]
finally_desire_cell <- what_want_from_cells[what_want_from_cells == "Significant"]

wanted_sected <- intersect(finally_desire_tis, finally_desire_cell)

length(wanted_sected) # 1
## [1] 1
venn_plot <- draw.pairwise.venn(length(finally_desire_cell), length(finally_desire_tis), length(wanted_sected), category = c("Tregs vs Tconv", "Breast tumor vs Normal breast tissue"), fill = c("blue", "pink"), ext.line.lty = "dashed")